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Abstract. We present results on the first excited states for the random-field 
Ising model. These are based on an exact algorithm, with which we study 
the excitation energies and the excitation sizes for two- and three-dimensional 
random- field Ising systems with a Gaussian distribution of the random fields. 
Our algorithm is based on an approach of Frontera and Vives which, in some 
cases, does not yield the true first excited states. Using the corrected algorithm, 
we find that the order-disorder phase transition for three dimensions is visible 
via crossings of the excitations-energy curves for different system sizes, while in 
two-dimensions these crossings converge to zero disorder. Furthermore, we obtain 
in three dimensions a fractal dimension of the excitations cluster of d s = 2.42(2). 
We also provide analytical droplet arguments to understand the behavior of the 
excitation energies for small and large disorder as well as close to the critical 
point. 
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1. Introduction 

Many disordered systems are characterized by an extremely rough energy landscape 
[H EJ [3] . This leads to the existence of configurations that are quite different from 
ground-state configurations but have an energy very close to the ground-state energy. 
The behavior of such systems differs considerably from the behavior of ordered systems 
and has attracted the attention of many physicists during the last decades. The most 
prominent example of a magnetic system with quenched disorder exhibiting many 
valleys in the energy landscape is the spin-glass model 0]. 

Another basic model from statistical physics that exhibits quenched disorder is 
the random-field Ising model (RFIM). In the RFIM, the sensitivity to the actual 
configuration of the disorder and the roughness of the energy landscape are not as 
strong as in spin glasses but nevertheless present and important to understand the 
behavior [5]. 

A RFIM sample consists of N = L d Ising spins s, on a d-dimensional regular 
hypercubic lattice with lateral length L in each dimension. Neighbouring spins interact 
ferromagnetically with strength J. In addition, on each site a random magnetic field 
hi acts on spin Sj. The hi are quenched random variables chosen according to some 
probability distribution. One common choice is a Gaussian normal distribution of zero 
mean and standard deviation h. With a Gaussian disorder, the ground state and the 
FES are nondegenerate. 

The Hamiltonian of the RFIM is 



The current state of knowledge on the Gaussian RFIM can be sketched as follows: 
No ordered phase can persist in the infinite system for d < 3, the lower critical 
dimension is di = 2 [6]. For d = 3, there is a continous phase transition at a critical 
field h c between an ordered ferromagnetic phase and a disordered paramagnetic phase 
[7] . It is possible to analyse the phase transition at zero temperature since the random 
fields are a relevant perturbation as shown by Fishman and Aharony [8]. 

Compared to other disordered models such as spin glasses, the RFIM has the 
practical advantage that in any dimension the problem of finding the ground state 
can be solved with fast graph-theory based algorithms. This means that the running 
time increases only like a low-order polynomial in N . This has led to numerous studies 
of the ground state landscape of the RFIM, e.g. [71 [HI US] • The system sizes that can 
be analyzed this way are considerably larger than those that can be equilibrated via 
Monte-Carlo techniques at small temperatures [TTJ fTJ] . 

Using these approaches one is not restricted to study only ground states. Even 
better, it is also possible to study the properties of the RFIM at very low but nonzero 
temperatures with suitably modified ground-state techniques. To do this, the ground 
state must be perturbed in some way to gain low energy states. Various attempts 
have been made in this direction, among them Refs. [3 [7]. 

In this study, we investigate first excited states (FESs) in the RFIM and their 
differences to ground states. In particular we present in Sec. [2] an algorithm to obtain 
exact FESs, in contrast to a previous attempt [13] . where some FESs were missed. 
In Sec. [31 we show first that indeed a considerable fraction of true FESs was missed 
in Ref. [13] ■ Next we apply our algorithm to study FESs both in the ordered and 
disordered phase. In particular, we investigate the scaling behavior of the excitation 
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energy and the fractal properties of the FESs. These findings are finally contrasted 
with arguments based on extreme-value considerations and on the droplet theory for 
the excitations in a hierarchical landscape. Finally, we finish with conclusions. 
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2. Algorithm 

A fast ground-state algorithm for the RFIM was first suggested in Ref. [14] and is 
explained in detail in Ref. [15] . The main idea is to transform the RFIM to a network 
Af = (G,c,s,t) consisting of a directed, weighted graph G = (V,E). V is the set 
of nodes and E C V x V being the set of edges. Cij > are edge labels denoting 
capacities. There are two distinguished nodes s and t called the source and the sink. 
All other nodes i ^ s,i ^ t are called inner nodes. In Af, a (s,t)-cut X = (S,S) is 
defined as a partition of V into two disjoint sets S and S with s G S and t G S. We 
call S and S the two sides of the cut. To each X, a capacity 

C(X)= ( 2 ) 
iesjes 

is assigned. Note, this means that only edges with the "right" direction i — ► j 
contribute to C(X). We will also use an equivalent description of a cut by a vector 
Xi G {0, 1} (i G V) with Xi = if i G S for X and X t = 1 else. 

The problem of finding the ground state of the RFIM can be reduced to the 
problem of finding the cut X mm — X mm (Af) with the smallest capacity among all 
possible cuts in such a network. This problem can be solved by determining the 
maximum flow of the network, as proven by Ford and Fulkerson 16 . As already 
mentioned, the ground state is nondegenerate, hence X mm is also unique for the 
model with Gaussian distribution of the random fields. 

We will now describe the algorithm of finding the first excited state as developed 
by Frontera and Vives (FV), adopting most of their terminology. By means of an 
example, we show why there are cases in which this algorithm fails to find the FES. 
The necessary adjustments are subsequently explained. 

Let Af be a network corresponding to a RFIM sample. Its minimum cut X mm 
can be calculated with a maximum-flow algorithm. We say an edge w — (i,j) 
is contained in A, if i G S and j G S. We denote the set of edges that are 
contained in X as E C (X) = G E\i G S,j G 5} and write w G E C (X). The 

set of nodes that are on either side of any edge w G E C (X) is called V C (X), i.e. 
V C (X) = {ie V\3j G V : G E C (X) or G E C (X)}. 

The set of all possible cuts (corresponding to all possible spin configurations) of 
Af is . Any edge w divides into two disjoint subsets, the set of cuts that 
contain w and the set S^f that do not. 

Next, a second network Af^ — (G, c%, s, t) is defined that includes a perturbation 
of strength S. The difference between Af& and Af is that the capacity of one edge 
w G E c (X mm ) of the minimum cut is increased by an amount 6 ^> 0, so that 

ci(w)=c(w) + 5 

c%{w) — c(w) yw / w. (3) 

Each cut X of Af is also well defined in Af^, but its capacity has increased if w is 
contained in X. For the remaining cuts of S£[ that do not contain w, the capacity 
does not change, i.e. 

CUX)-C(X) = { Hl X x (4) 

Since we chose w G E c (X mm ), X mm is the cut with the minimum capacity in 5^ 
for both A/" and A/"J (the capacity of all cuts of is increased by 8 in Af^). The 
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For X mm and all other 



(5) 



the maximum flow in the ground state), 



minimum cut of the complementary set is called X 
cuts of this set, the capacity is the same in both networks. 

It follows that depending on 6, there are two possibilities for the minimum cut of J\f^ 

if C{X min ) +6 < C(X 
if C(X min ) +6> C(X 

If S is chosen large enough (e.g. 5 = /^f, 
the minimum cut of A/"^ is always X mm . 

Let X 2nd e be the cut with the second smallest capacity in Af, corresponding 
to the first excited state. This cut must differ from the minimum cut by at least one 
edge. Frontera and Vives claim that because of this, "there exists (at least) one edge 
w of X mm so that X 2nd is in S^f" . We will show below that this assumption is not 
always true. 

To understand the algorithm of Frontera and Vives, we assume for a moment that 
the assumption is true. Then, in order to find the capacity of the second minimum 
cut, one would have to choose a suitably large 5. After that, C(X mm ) is determined 
for all networks Af^ depending on the edges w contained in the minimum cut X mm . 
From all the resulting cuts, the one with the minimum capacity in the original network 
M is accepted and called C*. 

C* = mh^ eEa(Xmm) C (X mm ) . (6) 

However, the FES has not necessarily been found yet, C* (j^nd smce ^he 
algorithm still has a problem: It assumes that the minimum cut necessarily includes 
an edge that the second minimum cut does not include. This is wrong, since there is 
also the possibility that all edges of the minimum cut are still contained in the second 
minimum cut which differs from the old one in that it additionally contains another 
edge. 

This can be illustrated by looking at Fig.Hl an example network where the GS cut 
is drawn in. When the capacity of the edge connecting s with node 1 is enlarged by a 
sufficiently large amount 5, the flow through the new network is limited by the edges to 
the sink. Therefore, in the resulting cut X* that is shown in the figure, all inner nodes 
have changed their side of the cut, and the capacity of X* is C* = 2 + 1.5 + 1.5 = 5. 
X* is not identical with the true second minimum cut C mm that is depicted in Fig. [2] 
and has a capacity C min = 3 + 1.5 = 4.5, so that C rnin < C*. 




Figure 1. Minimimum cut X mm and X* for an example network 
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Figure 2. Second minimum cut X 2nd for the example RFIM. 

We change the algorithm by leaving the original algorithm untouched (obtaining 
an intermediate cut C* after the first step) and introducing a second step. In this 
step, we go through all nodes k G V except those belonging to V c (X mm ). For each of 
these nodes k, we consecutively check if it is favorable to assign k to the other side of 
the cut (in the vector notation of a cut, set X^ ew = if X™ %n = 1 and vice versa). 
This would result in a new cut X new for which E c (X rnin ) C E c (X new ). As a result, 
additional edges contribute to the capacity of X new . We denote their contribution 
with C add so that C new = C min + C add . Note that performing this test for a node 
belonging to V c (X mm ) would result also in removing some edges from the cut, hence 
this is already taken care of by the FV algorithm. 

If there are nodes k where C new < C* , the FES is given by the minimum cut 
among the nodes k, if no such node exists, the FES is X* . 

Note that there can never be a negative contribution to C add since with > 
Vi,j this would mean to remove an edge that was included in X mm in the GS 
(while adding other edges at the same time) . But in this case the removed edge would 
have been already considered in the first step. 

For the extension of the FV algorithm, it is sufficient to look at single spin flips, 
i.e. at moving exactly one node from one side of the cut to the other. We explain this 
for a node k £ V with XJ! lm = 0, the other case X™ ln = 1 can be treated analogously. 
The change X™ in = -> X[! ew = 1 for k <£ V c (X mm ) affects the capacity of the new 
cut in the following way: 

There are possibly nodes X™ n = and Cki > 0. In this case an amount of cm is 
added to C add as soon as k is set to 1. Since k ^ V c (X mm ) there are no predecessor 
nodes i with Xf 1 ™ =1, hence incoming edges into k do not have to be considered for 
the calculation of C add , due to ©. 

It might seem plausible that there exist situations where, instead of adding the 
capacities cm to C add , also node i (and possibly more nodes) should be moved to the 
other side of the cut. In fact, some of these situations would lead to a smaller C add . 
However, in such a situation it would lead to an even smaller C add (remember that 
we are only interested in the second smallest cut) to invert only node i instead of 
inverting both node k and node i. This is due to the fact that with our method, there 
can never be a negative contribution to C add . 

The pseudo code for the full algorithm for calculating the FES is shown in Fig. 

El 
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algorithm First Excited State 

Transform the RFIM to a network Af 

Calculate X mm of Af that corresponds to the ground state 

C 2nd ._ (jmin + g with g = jGS^ 

for each edge w £ X mm { Frontera-Vives algorithm } 
c(w) := c(w) + S 
X w := minimum cut of Af^ 
if C'i(X w ) < C 2nd then 
C 2nd :=C S JX W ) 

v2nd v 

end if 

c(w) :— c(w) — S 
end for 

for all nodes k £ V c (X mm ) { extension to obtain true FES } 

£jadd q 

if X^ m = 1 then 

for all nodes i with Cki > 
C add ._ C add + Ch 

end for 

else if Xl lln = then 

for all nodes i with Cik > 
C add ._ C add + c , k 

end for 
end if 

|£ (jmin _j_ Qadd ^ £j2nd 
£j2nd ^jrain _|_ Qadd 

X 2nd ._ X min w j th node fc inverted 

end if 
end for 

return C 2nd and X 2nd 
end 

Figure 3. Pseudocode of the full algorithm to obtain the exact FES. 



Finally, we discuss the running time of the algorithm. The running time of the 
FV algorithm scales like [13] t ~ N a+1 when the running time of the minimum cut 
has scales as t ~ N a . The additional part introduced here is linear in N and therefore 
does not increase the full running time t of the algorithm in a noticeable way. 
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We calculated first excited states of 10 4 samples at different random field strengths h 
and different system sizes up to L = 16 in d = 3 (resp. L = 32 in d = 2). For the 
largest three-dimensional system size L = 20, 2 • 10 3 samples were calculated at each 
h. 

A comparison of the states generated with the Frontera-Vives algorithm with 
those generated using the complete algorithm that is described above, shows that 
the differences are noticeable though not drastic. For small systems in d = 2 with 
N = 4 3 spins at h = 2.3, the FV algorithm does not find the true FES in 15% of 
the cases. With growing system size, this percentage decreases. In two-dimensional 
32 x 32 samples (systems of this size were mainly investigated in [13]), the percentage 
is 5% and in three-dimensional systems of L = 16, 1% of the excitations generated by 
the FV algorithm are incorrect - both in the ferromagnetic and paramagnetic region. 
In the transition region between the two phases, the error rate is considerably smaller 
since a large number of first excited states consist of more than one spin and are 
therefore spotted correctly by the FV algorithm. 

Since the FV algorithm creates states with a higher excitation energy, the 
difference is also visible in the distributions of the excitation energy. The difference 
P{E) — P(E) fv is positive for small E and negative for larger E, as shown in Fig. [4] 

In our analysis of the properties of the FESs, the main quantities of interest are 
the energy difference E to the GS and the size V of the excitation cluster, given by the 
number of spins that change their orientation in the FES. By plotting the disorder- 
averaged excitation energy E against h for various L, we see that for large h, the E(h) 
becomes independent of h, i.e. E(h) — > const(i). The value of the plateau scales like 
E(h -> oo) ~ L~ d . 
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Figure 4. Difference P(E) — P{E)py between excitation energy distributions 
for the complete algorithm and the FV algorithm. Lines are guides to the eyes 
only. 
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d = 3 is shown in the main part of Fig. O the corresponding curve for d = 2 is plotted 
in the inset. 
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Figure 5. Excitation energy versus random field strength. Inset: E ■ L? versus h 
for d = 2. Main part: E ■ L 3 versus h for d = 3. Lines are guides to the eyes only. 

For h ^> h c , the curves of all system sizes coincide as expected. Slightly above 
the phase transition, the rescaled energy grows with decreasing h. For larger systems, 
the slope is larger. An interesting point is that in d = 3 the different curves intersect 
in a narrow region 2.1 < h < 2.3 that is close to h c . This behavior can only be seen in 
three dimensions. In d = 2 there is no single intersection point and the intersections 
approache h = with increasing L, reflecting the absence of a phase transition in 
two-dimensions [6]. Hence, it seems to be possible to locate the phase transition from 
the scaling behavior of the excitation energy. We are not aware of previous similar 
findings. In the next section, we will show by simple droplet arguments, that indeed 
the energy close to h c should display a L~ d behavior. Nevertheless, we are not able 
to give a full explanation of this behavior around h c , which seems to originate from 
corrections to scaling. 

The energy distributions P(E) are shown in Fig. [6] for fixed L = 16. Close to and 
beyond h c the distribution seems to be exponential, as we also find using the simple 
droplet arguments in the next section. In the inset, the extremely ferromagnetic case 
h = 1.0 is depicted. P(E) is peaked at a large E in completely ferromagnetic samples, 
corresponding to flipping spins against the orientations of all neighbors. For a detailed 
discussion, see again below. Small excitations cannot occur, there is an energy gap 
between the GS and the FES. 

The average size V (i.e. the number of flipped spins) of the FES is shown in Fig. 
[7] Global flips were not taken into account for this plot. In the ferromagnetic phase, 
V is equal to one. For h — *■ oo, it also decreases slowly towards V = 1. Thus, in both 
cases, the FES behaves as expected above. In between, there is a maximum with a 
size that seems to grow with L. The value of h where the maximum is located depends 
on L. It shifts towards h c w 2.27 with growing system size. This is due to the fact 
that finite systems are most "critical" at a field h c (L) > h c that approaches h c when 
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Figure 6. Probability distributions of the 3d excitation energy at fixed system 
size L = 16 for various h. In the inset, the range of the x-axis is [0:8], in the 
main part, the region of small E in the range [0 : 0.8] is shown. Lines are guides 
to the eyes only. 




Figure 7. Average volume of the FES dependent on h for various L. Lines are 
guides the eyes only. 



L — > oo. Since for first excited states, we are restricted to rather small system sizes 
and the statistical errors are still significant, it remains unclear if the maximum really 
moves to h c for larger systems. 

The probability distribution of the excitation volumes is shown in the double 
logarithmic plot of Fig. [H Both in very paramagnetic (h = 5.0) and ferromagnetic 
(h = 1.0) systems, P(V) declines in a fast way. At h = 3.0 near the peak of Fig. [71 
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Figure 8. Probability distribution P(V) at fixed L = 16 and random-field 
strengths h = 2.0, h = 3.0 and h = 5.0. Lines are guides to the eyes only. 

the decay is not so fast and resembles at least for small V a power law. However, to 
verify that P{V) approaches a power law at h c for large systems, much larger systems 
would be necessary. 

We also analyzed the shape of the excited clusters. The surface A of a cluster 
is given by the number of bonds that connect a spin of the cluster with another spin 
which is not part of the cluster. By plotting A against the volume V of the clusters, 
we gain information on the fractal dimension d s of the cluster. Since we found that 
the clusters are compact (V ~ R 3 , not shown) it follows that A ~ V ds ^ d . From Fig. 
M (h = 3.0, L = 12) it follows that d s = 0.808, leading to d s = 2.43(2) which can be 
compared with an earlier result d s = 2.30(4) of Middleton and Fisher [7] obtained for 
domain- walls. 

4. Extreme-value and droplet arguments for the minimal excitations 

Now, we try to understand the behavior of the distribution of the droplet energy and 
the average via analytical arguments for the extreme cases h <C /i c and h ^> h c as well 
as in the region where the droplet theory holds, hence close to h c . 

The relation E ~ iV™ 1 for large h can be explained by a simple extreme- value 
argument, which e.g. has also been used to explain the low-temperature behavior of 
spin glasses [17]. If h — * oo, the system behaves paramagnetically. Hence, FESs will 
be essentially single-spin flip excitations which is indeed true for the simulation data, 
as we have seen above. In this case, the spin with the smallest total local field will 
flip for the FES. Due to the paramagnetic state, there will be many spins with the 
same number of "up" neighbors and "down" neighbors. Hence, the total local field 
will be the same as the quenched local field hi. The energy of the FES will be given 
by the smallest absolute value \hi\ among these spinsU Thus, the excitation energy 

\ There might be also spins with the number of "up" spins not equal to the number of "down" spins, 
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Figure 9. Fractal dimension d s of first excitation droplet (L = 12) at h = 3.0. 
Each droplet is characterized by a data point. A least-square fit yields d s = 
2.43(2). 



/ r\h.tt\ \ 
P N (\h eff \) = NP(\h eff \) (l- J o dhP{h)\ 



is given by the minimum of O(N) absolute values of Gaussian random numbers. For 
simplicity, we assume O(N) = N: 

\h eff \ = \E\/2 = wm{\h 1 \ 1 \h2\ 1 ...,\h ff \} (7) 

The distribution function Pw(\h e ff |) of the sample minimum of N numbers from a 
distribution P(h) with P{h) = for h < is 

N-l 

(8) 

If P(h) is continuous and has a finite weight P(0) at the origin (which is true for a 
Gaussian distribution), one can approximate P(h) by P(0) for large N and gets 

p N (\h eff \) « ivp(o) (i - pco)!^//!)^ -1 

= NP(0) exp ((JV - 1) ln(l - P(0)\h eff \)) 
*\(N)exp{-\(N)\h eff \} 

with X(N) = P(0)N, having the mean E = T^(\h eff \) = i - TV -1 

P ~ A^ -1 agrees with the numerical simulation data of the mean energy and the 
predicted exponential decay of P(E) occurs indeed in the numerical data and can be 
seen in Fig. [6] 

The limit of strongly ferromagnetic systems can also be understood. In the 
ferromagnetic phase, there are two possibilities for the FES: 

In smaller samples, the entire system can flip. This happens if the excess of the random 
fields hi is even smaller than the minimum bond energy that has to be invested for 
flipping a single spin. If all spins flip, the bond contribution of the total energy stays 



but with the total local field close to zero. Nevertheless, since the Gaussian is peaked at zero, those 
will much less frequently comprise the FES compared to the spins discussed here. 
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invariant and the random field contribution E^ = J^. hi ~ L d / 2 changes its sign, so 
that the excitation energy is E = 2Eh- However, this type of excitation cannot persist 
in large systems, since Eh ~ £ 3 / 2 becomes soon larger than the energy of a second 
type of excitations, the local reversal of a single spin. This energy decreases for larger 
systems: 

If the system is fully magnetized, all spins are pointing in the same direction. 
We discuss the case that all spins are pointing "down" , the reverse case is similar. 
Hence, we have 2\hi\ < 12 for all spins. Then, from the y spins where the local field 
is positive, the spin with the maximum local field max{/ii} will have the smallest local 
field, i.e. will be chosen to be flipped for the FES. Since six bonds of strength J from Si 
to its neighbors are no longer satisfied, the excitation energy is E = 12 J — 2 max{/ii} 
in three dimensions. So we seek for the sample maximum in the ferromagnetic case. 

The sample maximum of N values drawn from a distribution P{hi) has the 
distribution 

)N-l 
(9) 

We calculated the integrals numerically and determined the expected value 
PN(hmax) for sample sizes up to N = 1000. It was found that PN(h max ) grows 
in a logarithmic way with N. Therefore we expect that 12 — E will also grow in a 
logarithmic way in systems that are large enough that global spin flips are negligible. 

The comparison with the numerical data of the FES deep in the ferromagnetic 
phase at h = 1.0 is shown in Fig. [TTJ] In the inset, the probability of a global flip is 
plotted versus N. This probability decreases quickly, and for L = 8, global flips occur 
already in less than 5% of the samples. In the main part, 12 — E is plotted versus 
N. 12 — E grows in a logarithmic way as expected, and deviations occur only for the 
smallest samples where global spin flips are still frequent. 



P N (h max ) = NP{K 
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Figure 10. Inset: Probability of global spin flips dependent on the number of 
spins N. Main part: 12 — E versus TV. The line represents a simple logarithmic 
function. 
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The droplet argument for the minimal excitation, which should be valid close 
to h c , starts with the assumption that the contributions from various scales are 
independent. Further, one takes their (scale-dependent) energy density to scale as 
a/L , where 9 is the energy fluctuation exponent [IS]. The value in 3D is currently 
found to be by ground-state simulations to be 9 = 1.49 ± 0.03 [7]. On scale I = 2 l , 
for small energies E, one will have thus the cumulative excitation energy distribution 
Pi = aE/l 9 . 

The smallest excitation energy in the whole system is then given by the mini- 
mum over all the I — 1 . . . L = 2 M scales. At each scale, the system splits into 
N(l) = (L/l) D independent subvolumes. The minimum excitation energy from a 
scale I has the cumulative distribution 

P(E,l) = l-(l-Pi) n W (10) 

where n(l) = (L/l) D . This comes simply from the argument that the complement of 
P(E, I) indicates that all the subvolumes have a smallest excitation of higher energy 
than E. Similar ways of reasoning have been used in other disordered systems [19] , 

The minimal energy (cumulative) distribution over all the scales can be 
constructed similarly. The droplet enery is smaller or equal to E if the converse 
is not true: that all the scales would have a minimum higher than that. Thus we 
obtain 

p(e) = i - p(e, o) = i - H(i - Pi) n(i) . (ii) 
i i 

The product is now taken over logarithmic scales, so that 1 < I < e lnL . 

The product can be converted by standard tricks to a sum and then to an integral, 
or computed numerically from 

P(E) = l-exp n(l)ln(l-Pi) (12) 

1=1... lnL/ln2 

where one notes that the argument of the exponential has a prefactor L d coming from 
the n(l). 

The resulting definite integral, using I — 2 X = e xln2 (x = 0, . . . , InL/ ln2), is 
of the type J exp (—ax) ln(l — frexp (— cx))dx and does not appear to have a simple 
closed form solution. Integrating, with fixed L over I = 1 .. .InL/ In 2 produces for 
the probability Figure [TTJ We see that in analogy to the numerical findings for 
paramagnetic 3d systems the minimum excitation energy distribution is exponential. 
Here, a in Pi is just chosen at random to a small constant. In reality, when the droplet 
argument is valid it should depend on dimensionality and h. Figure [12] shows for a 
range of L values the resulting average droplet energy E. In agreement with the most 
clear-cut numerical data in 3d, the scaling is volumelike (E ~ 1/L 3 ). 

5. Summary and discussion 

In this paper, we have first presented an algorithm that calculates first excited states 
in the random-field Ising model. We showed that a former algorithm developed by 
Frontera and Vives does not find the true first excited states in all cases but and 
therefore needs to be extended. When comparing the correct algorithm with the 
Frontera- Vies algorithm, we find that the difference in quantities such as the excitation 
energy is measurable albeit not drastic for larger systems. 
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Figure 11. Energy distributions for two fixed system sizes, computed via the 
assumption that one can equate L = 2 X where X is the number of independent 
lengthscales (corresponding to the 3d case). 




Figure 12. Average minimum droplet energy from the droplet argument, 
computed via integrating numerically Eq. I|12| l and the resulting average for the 
3d case. The data implies an inverse volume-like scaling E ~ 1/V = 1/L 3 . 

The excitations consist of connected clusters of spins which are reversed with 
respect to the ground state. In our simulations we study mainly the three-dimensional 
RFIM, but we also present some results for two dimensions for comparison. 
Interestingly, there is an intersection of the E(h, L)L 3 curves for different system 
sizes at a point very close to the phase transition point h c . This behavior cannot 
be found in the two-dimensional case and seems to reflect the existence of a phase 
transition. 

The volume distribution P(h) exhibits a peak close to h c that becomes sharper 
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with growing system size. However, the size distribution P(V) at a fixed h near the 
peak is relatively broad so that statistical errors are large. The fractal dimension of 
the clusters was determined to be d s — 2.42(3). 

Using analytical arguments, we were able to understand the behavior of the 
average excitation energy E(h) in the extreme cases h ^> c and h <C h c . To further 
investigate the origins of the observed energies and their scaling, we have used a droplet 
argument to derive the energy distribution of the minimal excitations. This line of 
reasoning is similar to a few others found in the literature (see again e.g. [15] . but 
we are not aware of any where the whole distribution would have been written down. 
The result reproduces well some features of the numerics, as a volume-like scaling of 
E and the exponential character of the P(E). Hence, the such calculated behavior 
at h c follows our numerical observations. Additional corrections to scaling close to 
h c seem to be responsible for the observed behavior of the crossing of the E(h, L)L 3 
curves. 

For future work, it would be interesting to look at other types of low excitations 
with ground states that do not necessarily produce the FES but still a low energy 
state. States of this kind can be generated much more efficiently, leading to larger 
systems that can be simulated. Another possibility would be to analyze excitations 
in higher dimensional systems (see [20] [H] for results in d — 4). However, different 
approaches from the FES algorithm would have to be used for that since we were 
already restricted to relatively small system sizes in three dimensions. 
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